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In this work, we study the in-vitro dynamics of the most malignant form of the primary brain 
tumor: Glioblastoma Multiforme. Typically, the growing tumor consists of the inner dense pro- 
liferating zone and the outer less dense invasive region. Experiments with different types of cells 
show qualitatively different behavior. Wild-type cells invade a spherically symmetric manner, but 
mutant cells are organized in tenuous branches. We formulate a model for this sort of growth using 
two coupled reaction-diffusion equations for the cell and nutrient concentrations. When the ratio of 
the nutrient and cell diffusion coefficients exceeds some critical value, the plane propagating front 
becomes unstable with respect to transversal perturbations. The instability threshold and the full 
phase-plane diagram in the parameter space are determined. The results are in a good agreement 
with experimental findings for the two types of cells. 



PACS numbers: 87.18.Ed, 87.18.Hf 

One of the most aggressive forms of primary brain tu- 
mor is Glioblastoma Multiforme (GBM) Despite ma- 
jor advances in medical science the prognosis for victims 
of this disease is very poor 0: the median survival for 
patients with newly diagnosed GBM is approximately 12 
months. One of the main reasons for such high mortal- 
ity and poor success of treatment is the fact that GBMs 
are highly invasive |2| . The growing tumor sheds invasive 
cells which run through the brain, see Fig.^ The invasive 
nature of malignant gliomas makes treatment difficult Q ; 
secondary tumors are produced by the invasive cells even 
if the primary is removed. In this paper we introduce a 
reaction-diffusion model for invasion. By comparing two 
different cell lines we hope to get insight into invasion 
dynamics which is needs to be better understood. 

This work is inspired by recent in vitro experiments 
0, |(| where microscopic tumor spheroids (radius about 
250 n) were placed in collagen-I gel and allowed to grow. 
The cell lines used were U87 and U87-AEGFR. The first 
type is called 'wild- type' in what follows. The second 
is a mutant line 4] in which there is an amplification of 
the epidermal growth factor receptor (EGFR) gene. This 
amplification occurs in approximately 40 percent of cases 
of GBM @. 

If we compare the growth of the two cell lines in vitro 
we see two main differences; cf. Figure ^ The inva- 
sive region for the wild-type cells grows faster than for 
the mutant cells, and the wild-type produces a spheri- 
cally symmetric pattern, whereas mutant cells produce a 
branching pattern [(| . In the present work, we formulate 
a simple reaction-diffusion model that is able to repro- 
duce these experimental findings and may give insight 
into the functional significance of the mutation. 

We interpret the branching shown in Fig. [Tfe, as a 
branching instability. Analogous instabilities were stud- 
ied in the theory of combustion |7| , and in studies of the 
self-organization of microorganisms || • One way of mod- 
cling this is to assume that there is attraction between 




FIG. 1: Growing tumors from in vitro experiments |y in col- 
lagen gel for the wild-type (a) and mutant (b) cells. These 
in vitro tumors consist of an inner proliferation zone with a 
very high density of cells and an outer invasive zone, where 
the cell density is smaller. The structure in vivo is believed 
to be similar. The radius of the inner zone here is about 250 
u. For wild-type cells a spherically symmetric pattern is ob- 
served, (a). Mutant cells are organized in tenuous branches, 
(b). Note also that the invasive region for the mutant type 
cells grows slower than for wild-type cells. 



cells [9j due to the production of growth factors. Another 
possibility is to assume nonlinear diffusion, where the dif- 
fusion coefficient increases with the density of the cells, 
as was proposed for bacterial colonies We will intro- 
duce another mechanism based on the known biology of 
GBM. 

In our continuum description, we will deal with the 
density of cells u(r, t) and the density of some growth 
factor or nutrient (whichever controls the growth), c(r, t). 
We assume that c diffuses to the tumor from far away. 
Each cancer cell is able to proliferate as well as to per- 
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form random motion. It is known that within the inner 
region, cells have quite a high proliferation rate whereas 
in the invasive re gion , cells have high motility, but low 
proliferation rate |lOj . This is an indication of a dynam- 
ical switch between the cell phenotypes 0. We model 
it by introducing a density- dependent proliferation term 
where the proliferation rate increases with cell density. 
The simplest form for such a term is du/dt oc u 2 c. (The 
extra power of u compared to the usual uc means that 
high density gives rapid proliferation since the prolifera- 
tion rate per cell, (l/u)du/dt oc u.) As we will see, this 
term drives the instability and leads to branching. We 
assume that in order to proliferate a cell needs to con- 
sume some amount of c. The density c obeys a diffusion 
equation with a sink at the tumor cells. 

We encode these assumptions in the following equa- 
tions: 



du 

dt 
dc 

di : 



V- (D u Vu) + au 2 c, 
V- (D c Vc) - (3u 2 c. 



(1) 



Here D u and D c are the diffusion coefficients of u and 
c, a is a proliferation coefficient, and (3 is the coefficient 
of nutrient consumption. We assume that the density in 
the center of the tumor is not very high compared to the 
density of closely packed cells u c . We suppose also that 
the nutrient concentration is kept constant far from the 
tumor: lim r _, oo c(r) = c^. 

In what follows we will measure cell density in the 
units of some characteristic density uq, nutrient density 
in units of Coo, distance in units of [I? c /(/3tio 2 )] 1 ^ 2 , and 
time in units of (Puq 2 )^ 1 . This gives: 



— = Iv 2 

dt 5 

dc _ 2 
dt = VC 
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where S — D c /D u is the ratio of the nutrient and cell 
diffusion coefficients and m = (3uQ/{ac ao ) is the ratio of 
consumption and proliferation rates. 

Typically, the nutrient or growth factor represented 
by c is a small molecule. It is expected to diffuse much 
faster than the cells. For example, the diffusion coeffi- 
cient of glucose in the brain is of the order of 10~ 7 cm 2 /s, 
while the cell diffusion is of the order of 10~ 9 cm 2 /s |9|, 
so that 8 ~ 100. A typical nutrient consumption is 
10~ 12 g/cell/min 0], and a typical glucose concentration 
is of the order of lg/1. Assuming that typical cell density 
within the invasive region is of the order of 10 5 cell/cm 3 , 
we estimate the consumption rate as 1.7 x 10~ 6 S _1 . The 
typical proliferation rate in experiments is of the order 
of 1/day, so that m turns out to be of the order of 0.1. 

We work in a two-dimensional channel geometry. Let 
x be the direction of tumor growth, and y be the trans- 
verse direction, perpendicular to the direction of the front 



propagation. In the y direction we use periodic boundary 
conditions with a finite channel width. In the x direction, 
far ahead of the tumor, the cell concentration is zero, 
u(x = oo) = 0, and the scaled nutrient concentration is 
unity, c(x = oo) = 1. On the other hand, at x = — oo we 
demand c = 0. There is a conservation law in Eqs. J2J: 
a volume integral over the system of (mu + c) is a con- 
served quantity. Its interpretation is that in our model a 
cell needs some amount of food to divide. Therefore, at 
x = — oo, u = l/m if there is a steady state. Our initial 
conditions are u = l/m for x < 0; u = for x > 0, and 
c = 0, for x < 0; c = Coo for x > 0. We will investi- 
gate the situation after transients have died away and a 
steady propagating state has been established. 

First, we consider the solutions of Eqs. J2J in the form 
of plane propagating fronts: u — uq(^);c = co(£), £ = 
x — vt. Substituting into Eqs. (0 we arrive at: 
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To obtain the profiles and the velocity of front propa- 
gation, we performed (in Matlab) the following shooting 
procedure. First, we write down Eqs. © as four cou- 
pled first-order differential equations in the form (o)' = 
M(a), where a is the column of solutions with elements 
Mo, uo', Co, Co'. Then we find the eigenvalues and eigen- 
vectors of M at £ = ±oo. Starting with the solution a 
that is proportional to the eigenvector belonging to the 
positive eigenvalue at £ = — oo (using the linearity of the 
problem, we choose the constant of proportionality to be 
unity), we perform shooting by the velocity of front prop- 
agation v. We find the profiles by demanding that the 
solution a at £ = +oo is a linear combination of the two 
eigenvectors ipi and tp2 belonging to negative eigenvalues 
Ai and A2 at £ = +00: a = hy^\ exp(A 1 £)+& 2 V'2 exp(A 2 £). 
Figure [21 shows a typical solution of Eqs. J5J for u = 
and c = c(£). 
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FIG. 2: Density profiles of cells and nutrient from Eqs. 
The first curve is the cell density u = w(£) (solid line.) The 
second curve (dashed) is c = c(£) (inset) . The parameters are 
m = l,S= 100. 
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We now perform a linear stability analysis. Consider 
perturbations in the transverse direction 

u = Uo(£) +u 1 (£)exp(jt + iky), 

c = c (£) + a (0 exp ( 7 t + ifcy) (4) 

and substitute into Eqs. @. We have 

1 / 2 fc 2 \ 1 

-m" 7 vui + I — c u - -y - 7 u i H «o 2 ci = , 

o \m o J m 

c\" 7 wci' - (m 2 + fc 2 + 7) ci - 2c u ui = . (5) 

For a fixed value of transverse wave number k, we should 
find the perturbations and Ci(£), and the growth 

rate 7. As before, we rewrite Eqs. as four coupled 
first-order differential equations in the form (au n )' — 
Miinifliin). There are two positive eigenvalues of the 
matrix Mu n at £ = —00. We start at £ = —00 from 
the linear combination of the corresponding eigenvec- 
tors, aim = b^ips exp(A3£)+&4'!/'4 exp(A4^), where one can 
chose &4 to be unity due to the linearity of the problem. 
Performing shooting in two parameters, the growth rate 
7 and the constant 63, we find the eigenfunctions 
and Ci(£) by demanding that the solution 3u n is given 
at £ = 7 00 by a linear combination of the two eigenvec- 
tors ^5 and ipe belonging to negative eigenvalues A5 and 
A 6 , a Hn = b 5 ip 5 exp(A 5 £) + b 6 tp 6 exp(A 6 £). Changing the 
value of transverse wave number fc, for the fixed S and 
m, we calculate the dispersion curve 7(fc). 



** 










/4 


k, 




* 




/ 3 




* 
















2 




is-'' 






1 


/ 






V 




4 




e 




<«- 










0.5 


1 





FIG. 3: An example of the dispersion curve 7(fc). The in- 
stability occurs if the ratio of diffusion coefficients 8 exceeds 
a certain critical value. For S > 8 cr ~ 2.300, the growth 
rate 7 is positive for small k, while for larger fc, cell dif- 
fusion in the transverse direction stabilizes the instability. 
The parameters are: m = 0.1, S = 20. The inset shows 
the dependence of the largest unstable wave number fc* on 
e = (5 — 5 cr )/8 cr . Numerical simulations are the asterisks, 
the asymptote fc* = (0.36/m)e 1,/2 is the dotted line. 

As was found previously in the context of chemical re- 
actions for m = 1 |l2l | , plane fronts can become transver- 
sally unstable if the ratio of diffusion coefficients 5 ex- 
ceeds a certain critical value. Indeed, for 6 > 8 cr , the 
growth rate 7 is positive for small fc, while for larger 



fc, cell diffusion in the transverse direction stabilizes the 
instability. We checked that 5 cr « 2.300, in agreement 
with previous results 0,0. Figure shows an example 
of the dispersion curve for S > 6 cr . An inset shows the 
dependence of the largest unstable wave number fc* on 
e = (5 — S cr )/S cr . Numerical simulations are denoted by 
asterisks, the asymptote fc* = (0.36/m)e 1 / 2 is shown by 
the dashed line. 

For a very wide system the instability threshold S cr 
does not depend on m. To see this, introduce new di- 
mensionless variables r = {m/ 5 1 / 2 )R, t — m 2 T, and 
u = U/m. In this case, m drops out of the problem. 
A consequence of the elimination of m is that one can 
easily find the dependence of the velocity, v, and of the 
wave number, fc, on m: v = m~ 1 5~ 1 / 2 V, k = m~ 1 8 1 / 2 K. 
Since the scaled front velocity, V, and the scaled wave 
number, K, must be independent of m, v and fc are pro- 
portional to mT x . We will not eliminate the parameter 
m from the problem and will work with Eqs. |2"| l. Differ- 
ent types of cells have different diffusion coefficients D u 
and different proliferation rates a. Therefore, it is con- 
venient that the dependence of the physical quantities 
kphys = k[D c /(/3u 2 )}~ 1/2 and v phys = -y( J D c /3u 2 ) 1/2 on 
D u and a enters only via the dimensionless variables v 
and fc. 

For a system of finite width we need to form a discrete 
set of modes from the modes of the infinite- width system 
by imposing kL — 2mr, where n = 1,2,... and L is the 
width of the system. For a small enough system we can 
'freeze out' the instability if fc* < 2ir/L. For a spherical 
tumor L is of the order of the diameter. For example, 
Figure ^a), could correspond to small fc*. 

We now consider a phase plane of parameters (5, m), 
see Fig. 0] As mentioned above, the main differences 
between the experiments with wild-type and mutant cells 
are in the velocity of the front propagation and possible 
symmetry-breaking. It means that the region of larger v 
and smaller fc* in this phase plane corresponds to wild- 
type cells, while the region of smaller v and larger fc* 
corresponds to mutant cells, see Fig. 0] 

We consider two families of curves: the first is of con- 
stant front velocity v = const and the second is of con- 
stant largest unstable wave numbers fc* = const. We 
focus first on v = const curves. It was shown previously, 
that for m = 1 and large 5,v = 1.219<5 -1 0|. Combining 
this with the m dependence, one can see that in the (<5, m) 
phase plane the curves of constant velocities for large S 
are given by 5 = const/m. Then, we consider fc* = const 
curves. Our numerical calculations indicate that for large 
values of 6, the largest unstable wave number fc* tends 
to some constant independent of 6; the value of this con- 
stant for to = 1 is approximately 0.5. Therefore, for large 
5, fc* = 0. 5/m, and the curves of constant fc* are given 
by to = const. Figure 01 shows also these large S asymp- 
totes. A typical length scale is [D c / {f3uo 2 )] 1 / 2 ~ 0.2cm 
and a typical velocity scale (DcPuq 2 ) 1 ^ 2 ~4x 10~ 7 cm/s. 
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FIG. 4: The phase plane (<5, m) with two families of curves: 
v — const and fc» = const, calculated from Eqs. @ and 
The two curves marked by o's show v = const; v — 0.1 (left), 
and ii = 0.02 (right) . The two curves marked by *'s show 
fc, = const . The left-hand curve corresponds to a stronger 
instability, fc* = 10, while the right one corresponds to a 
weaker instability, fc* = 1. Also shown are the large 8 asymp- 
totes: 5 = const/m for v = const curves, and m = const for 
k- t — const curves, see text. The instability threshold for an 
infinite stripe, 5 ~ 2.300, is plotted by the dotted line. The 
larger-5 and smaller-m region correspond to wild-type cells, 
while mutant cells are in the region of smaller <5 and larger m. 
This suggests that wild-type cells have a larger diffusion con- 
stant, but a smaller proliferation rate, compared to mutant 
cells. 



Comparing this with experimental data [(|, one can see 
that the dimensionless wave number and front velocity 
v should be of the order of 5 and 0.1, correspondingly. 

As one can see on Fig. larger-<5 and smaller-m re- 
gion corresponds to wild-type cells, while mutant type 
cells correspond to smaller-<5 and larger-m region. Thus 
we predict that wild-type cells have a larger diffusion con- 
stant but a smaller proliferation rate than mutant cells 
in the invasive zone. Comparing these theoretical pre- 
dictions with experiments |6j , twopoints should be taken 
into account. First, experiments [6| with tumor cells were 
performed in a three-dimensional geometry, so the tumor 
growth can be described by spherical propagating fronts 
for wild-type cells. In order to explain branching patterns 
of mutant type cells, linear stability analysis of spherical 
fronts, rather than plane fronts, should be performed. In 
this case, the role of the initial tumor radius should be an- 
alyzed. Second, the basic solutions of our model are plane 
fronts which propagate with constant velocity. However, 
in experiments, the more dense inner proliferative region 
grows slower than the less dense outer invasive region Q . 
Probably the tumors are in a transient regime, and the 
steady-state behavior will set in later ^{|. However, our 
qualitative predictions should hold in three dimensions. 
Experimental verification of these features of the growth 
should be possible. 

In summary, we have considered the growth of GBM 
tumors. In vitro experiments 6] showed that the dynam- 



ics of growth and resulting patterns are quite different for 
wild-type and mutant cells. For the wild-type the inva- 
sive region grows faster, and tumor remains spherically 
symmetric. On the other hand, the invasive region grows 
slower for the mutant cells, and there are indications of 
symmetry-breaking of spherically symmetric growth. We 
formulated a simple reaction diffusion model that cap- 
tures these experimental findings. Based on our model, 
we explain different patterns by different diffusion con- 
stants and proliferation rates of wild-type and mutant 
cells: wild-type cells diffuse faster, but have a lower pro- 
liferation rate in the invasive zone. We think that an 
attempt should be made to test these predictions and re- 
late them to the microscopic biology of the two cell lines. 
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